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Abstract 

An exactly solvable phase-field theory of dislocation dynamics, strain hardening and hysteresis in 
ductile single crystals is developed. The theory accounts for: an arbitrary number and arrangement of 
dislocation lines over a slip plane; the long-range elastic interactions between dislocation lines; the core 
structure of the dislocations resulting from a piecewise quadratic Peierls potential; the interaction be- 
tween the dislocations and an applied resolved shear stress field; and the irreversible interactions with 
short-range obstacles and lattice friction, resulting in hardening, path dependency and hysteresis. A chief 
advantage of the present theory is that it is analytically tractable, in the sense that the complexity of the 
calculations may be reduced, with the aid of closed form analytical solutions, to the determination of the 
value of the phase field at point-obstacle sites. In particular, no numerical grid is required in calculations. 
The phase-field representation enables complex geometrical and topological transitions in the disloca- 
tion ensemble, including dislocation loop nucleation, bow-out, pinching, and the formation of Orowan 
loops. The theory also permits the consideration of obstacles of varying strengths and dislocation line- 
energy anisotropy. The theory predicts a range of behaviors which are in qualitative agreement with 
observation, including: hardening and dislocation multiplication in single slip under monotonic loading; 
the Bauschinger effect under reverse loading; the fading memory effect, whereby reverse yielding grad- 
ually eliminates the influence of previous loading; the evolution of the dislocation density under cycling 
loading, leading to characteristic 'butterfly' curves; and others. 
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1 Introduction 

This paper is concerned with the formulation of an exactly solvable phase-field theory of dislocation dy- 
namics, strain hardening and hysteresis in ductile single crystals. The present theory is an outgrowth of 
the phase-field model of crystallographic slip of Ortiz ([[]), which was restricted to a single dislocation loop 
moving through point obstacles of uniform strength and was based on a line-tension approximation. The 
present theory accounts for: an arbitrary number and arrangement of dislocation lines over a slip plane; the 
long-range elastic interactions between dislocation lines; the core structure of the dislocations described by 
means of a Peierls potential; the interaction between the dislocations and an applied resolved shear stress 
field; and the irreversible interactions with short-range obstacles and lattice friction, resulting in hardening, 
path dependency and hysteresis. 

A chief advantage of the present theory is that it is analytically tractable, in the sense that the complexity 
of the calculations may be reduced, with the aid of closed form analytical solutions, to the determination 
of the value of the phase field at point-obstacle sites. In particular, no numerical grid is required in cal- 
culations. On the sole basis of the value of the phase field at the obstacle sites, the theory reconstructs, 
analytically and in closed form, the location of — possibly large numbers of — discrete dislocation lines, of- 
ten in complex arrangements and undergoing intricate topological transitions, and endows each dislocation 
line with a well-defined core. Furthermore, the theory characterizes the equilibrium configurations of the 
dislocation ensemble directly, with the result that no short transients associated with artificial kinetics need 
to be resolved in calculations. By virtue of these attributes, the theory would appear to be advantageous in 
comparison to other recent phase-field models of dislocation dynamics (|2]; ||), which require the use of large 
computational grids and time integration over large numbers of time steps. 

For definiteness, in this paper we restrict our attention to systems consisting of a dislocation ensemble 
moving within a single slip plane through a random array of discrete forest dislocations under the action 
of an applied shear stress. In the present theory, the dislocation ensemble populating the slip plane is 
represented by means of a scalar phase field. Specifically, the value of the phase field at a point of the slip 
plane simply records the extent of slip in quanta of Burgers vector. The phase field is, therefore, integer 
valued, and its value at a point may alternatively be regarded as recording the number of dislocations, with 
proper accounting for their sign, which have crossed the point. Thus, the individual dislocation lines may 
be identified with the lines over which the phase field jumps by one. The phase-field approach furnishes 
a simple yet efficient means of representing arbitrary dislocation geometries, possibly involving very large 
numbers of individual dislocations, as well as enabling the tracking of intricate topological transitions, 
including loop nucleation, pinching and the formation of Orowan loops. 

The central objective of the theory is to characterize the evolution of the phase field. The theory is 
couched within a general variational framework for dissipative systems developed by Ortiz et al. ([|; ||; |6|; 
[7|; ||), and accounts for energetic and kinetic effects. The energy terms contemplated by the theory include: 
the core energy of the dislocations, represented by a piecewise-quadratic Peierls potential (§); the long- 
range elastic interactions between primary dislocations and between primary and forest dislocations; and 
the energy of interaction with the applied resolved shear stress. The particular piecewise quadratic form 
of the Peierls potential adopted here lends itself to an effective analytical treatment based on the Fourier 
transform @. The main effect of the Peierls potential is to endow the dislocations with a well-defined core 
energy and to coarse-grain the slip plane by suppressing wavelengths shorter than the lattice parameter. 

At zero temperature and in the absence of irreversible processes, the equilibrium configurations of 
the dislocation ensemble follow directly from energy minimization. The resulting variational problem is 
strongly nonlinear owing to the all-important constraint that the phase field be an integer-valued function. 
Furthermore, the minimization problem is nonconvex owing to the multiwell structure of the Peierls poten- 
tial. In addition, the variational problem is nonlocal, owing to the presence of long-range elastic interactions. 
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These attributes render the energy minimization problem mathematically non-trivial and confer its solutions 
a rich structure. Despite these difficulties, the choice of a piecewise quadratic Peierls potential lends the 
problem analytical tractability and, remarkably, a general solution to the problem can be found analytically 
in closed form. 

In the present theory, hysteresis arises from the assumed irreversible short-range interactions between 
primary and forest dislocations and from lattice friction. The strength of some of these interactions has 



recently been investigated using atomistic and continuum models ( JL0| ; |1 1| ; |12| ). Thus, we assume that the 
crossing of a forest-dislocation site by a primary dislocation always costs energy, regardless of the direction 
of the crossing. In the context of a phase-field representation, this is tantamount to assuming that any 
variation in the phase field at sites occupied by obstacles requires the supply of a certain amount of work, 
regardless of the sign of the variation. Physically, this work may be identified with the energy required 
to dissolve the dislocation-obstacle reaction product. Evidently, this form of interaction is irreversible 
and dissipative, and thus cannot be described by means of an energy function. It is possible, however, to 
develop an incremental variational framework which characterizes the evolution of the system by means of a 
sequence of minimization problems (Q; ||; ^; [7| |8|). The function to be minimized over a given step includes 
the incremental work of dissipation incurred as a result of obstacle crossings, and depends on the state of 
the system at the beginning of the step, which results in path-dependent and hysteretic behavior. 

The theory predicts a range of behaviors which are in qualitative agreement with observation, including: 
hardening and dislocation multiplication in single slip under monotonic loading; the Bauschinger effect 
under reverse loading; the fading memory effect, whereby reverse yielding gradually eliminates the influence 
of previous loading. Simultaneously with the deformation and hardening characteristics of the system, the 
theory naturally predicts the evolution of the dislocation density. In particular, no independent equation of 
evolution for the dislocation density needs to be supplied. 



2 Dislocation energies 

In this section we derive explicit expressions for the energy of a dislocation ensemble contained within a slip 
plane in an elastic crystal. The dislocation ensemble is described by means of a phase field defined over the 
slip plane. In addition to the elasticity of the crystal, the slip plane is endowed with a piecewise-quadratic 
Peierls interplanar potential. The elastic interaction energy formulated in this section is a special case of a 
general class of energies for continuously distributed dislocation loops in isotropic elastic crystals derived 



by Ortiz and Xu (13). Extensions to anisotropic crystals have also been given by Xu (14) 



We consider a crystal undergoing deformations characterized by a displacement field u,,. Following 



Kroner (15), we begin by decomposing the displacement gradient in the additive form: 



where uy is the displacement gradient, or distortion field, and, here and subsequently, commas are used to 
denote partial differentiation. For simplicity, we shall specifically focus on crystallographic slip occurring 
on a single slip plane S. Under these conditions, the plastic distortion j3fj is supported on S and has the 
form 

Pf j = S i m j 5 s (2) 

where Si is the displacement jump across S, m is the unit normal to S, and 5s is the Dirac distribution 
supported on S. The elastic distortion f3fj is assumed to be continuous. Neither the elastic nor the plastic 
distortion is required to be compatible, i. e., to be a gradient. 
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We shall additionally assume that the energy of the dislocation ensemble may be written as the sum 
of three terms: a core energy expressible in terms of a Peierls interplanar potential; the elastic interaction 
energy of the dislocations; and the energy of interaction with the applied stress field. These assumptions 
lead to the consideration of an energy functional of the form: 

E = j <P(S)dS + J -djuftjfarfx - J UdidS = E coie + E mt + E ext (3) 

Here, <fi denotes the Peierls interplanar potential and the first term in (||) represents the misfit or core energy 
of the dislocations. In the second term, which represents the elastic interaction energy, Cijki are the elastic 
moduli of the crystal. Finally, the third term represents the interaction energy between the dislocations and 
a self-equilibrated applied stress field resulting in a distribution of tractions t, over the slip plane. The aim 
now is to derive a phase-field representation of each of the terms in (||). 



2.1 Elastic interaction energy 

Suppose that the slip distribution Si over S, or equivalently, is prescribed. The corresponding elastic 
interaction energy may then be computed as follows. Insertion of ([]]) into (|3j> and minimization with respect 
to the displacement field yields the equilibrium equations 

(cijkiUk,i),j ~(cijkiPli),j = (4) 
The corresponding displacement field is: 

Uk — Gki ~k {pijmnf^rnn) ij (^) 

where Gki is the Green's function and (*) denotes the convolution operator. The distortion field then follows 

as 

Uk,l = —Gki,l * ( c ijmn(3mn)ij (6) 

and the elastic distortion as 

Pkl = ~Gki,l * ( c ijmnf3mn)ij ~Pkl 

Finally, the elastic interaction energy E mt is obtained by inserting <^> into the second term of @. 

A convenient explicit expression for E mt may be obtained by recourse to the Fourier transform. Using 
the convolution theorem, the Fourier transform of the elastic distortion follows from (^) as 

Phi = G kikjkiCij mn 0^ nn — (3^i (8) 

where a superposed Q denotes the Fourier transform of a function, hi is the wavenumber vector, and the 
Fourier transform of the Green's function is determined by the relation: 

G ik = c ijkikjki (9) 
Finally, the interaction energy follows from an application of Parseval's identity, with the result: 
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Here the symbol (*) denotes complex conjugation and we write 

-Amnuvfe) Cklr sCijmnCpquv(G jk-l ^klij^ r P k Q ks ^rspq) (H) 

— Cmnuv CkluvCijmnGkikjkl (12) 

Mura ( |l6| ) has shown that E mt can also be expressed in terms of Nye's (17) dislocation density tensor. 
Indeed, it is readily verified from <^> and @ that, if the plastic distortion [5f- is compatible, i. e., if it is a 
gradient, then the elastic distortions, and hence the elastic interaction energy, vanish identically. 

In order to facilitate analysis we shall resort to several simplifying assumptions. Thus, for the special 
case of an isotropic crystal the Green's function reduces to 

a- = i- ("a. - ( i3) 



13 2ix\k 2 1 - v A; 4 j 

where \i is the shear modulus of the crystal, v is its Poisson's ratio, and we write k = \k\. The interaction 
energy (10) then simplifies to 

Eint= jtf / 2 {^-(v-rnndf-iS^f + ^iS-^^rj.rnf^d'k (14) 

where we write r\i = ki/k. In addition, we shall adopt the constrained displacement hypothesis of Rice 
(|j"8|; p^ ) according to which displacements, and the attendant shear resistance, take place predominantly in 
the direction of the dominant Burgers vector, i. e., 

6i(x) = 5(x)si (15) 

where s = b/\b\ is a constant unit vector and 5 is a scalar function. For certain crystals, this conjecture has 
found support in atomistic calculations (20; 21 ; 22; 23). By virtue of these assumptions, ( fl4] ) simplifies to 

E mt = 7^-J^U + -Arvht) \S\ 2 d 3 k (16) 



(2vr) 3 J 2 

where we have chosen axes such that the slip plane is parallel to the (x\, X2)-plane and the Burgers vector 
points in the direction of the j;i-axis. Finally, for the particular case in which the slip distribution is confined 
to the plane X3 = 0, the elastic interaction energy further reduces to 

R int = 1 [ ft I k 2 , _J k l \ \£\2 H 2 k (l7) 

(2*)* J 4 l vv ^Tfcf + i-- VWTMJ 11 () 

which is the sought expression. 

By way of a simple illustrative and verification example we may consider the case of an edge dislocation 
pile-up, for which the slip distribution is of the form 

5 = 27rb6(k 2 )f(k 1 ) (18) 



whence the elastic energy per unit length of dislocation follows from ( |17| ) as 

E int fib 2 



L 4vr(l - v) 



hWfih^dh (19) 



5 



M. Koslowski, A. M. Cuitiho and M. Ortiz 



An Exactly Solvable Phase-Field Theory. 



For a single Volterra dislocation at the origin, f(kx) = \jik\ and the elastic energy per unit length evaluates 
to 

E int _ fib 2 W r ° dki __ fib 2 R 

~T " 4tt(1 - u) J„ /R "fcT " 4tt(1 - v) ° g ~ U)) 

where, in order to avoid logarithmic divergences, we have introduced lower and upper cutoff radii. We 
verify that (EO) coincides with the well-known expression for the energy per unit length of a straight edge 



dislocation in an isotropic crystal (24). For a screw pile-up, a similar derivation yields: 

pint ,,i,2 >>oo 



E mt fib 1 



\k 2 \\f(k 2 )\ 2 dk 2 (21) 



L 4tt 

For a single screw dislocation at the origin the dislocation energy per unit length evaluates to 

gmt «r/ro ^ ^2 R 

L 4vr J n/R k 2 4vr r 



which again coincides with the classical result (24). 



2.2 Core energy 

In the Peierls theory of the dislocation core, the interplanar potential (f>(5) is identified with the energy per 
unit area that results when two semi-infinite crystals are taken through a relative rigid displacement Si. From 
symmetry considerations it follows that, in the absence of an applied field, the energy of the crystal attains 
minima when the displacement jump Si is an integral multiple of a Burgers vector b{ of the lattice, i. e., at 

<5i = £&i, £GZ (23) 

where Z denotes the set of all integer numbers. These special slips determine the location of the wells of (f). 
With a view to enabling the application of Fourier transform methods, we shall assume that <p is piecewise 
quadratic (Q). In this model, the interplanar potential is taken to be of the form 

0(5) = mm -CjitSj - ZbjXSt - ih) (24) 



The moduli Cji may be determined by equating (g4|) with the energy per interatomic plane of a crystal 
undergoing a simple shear deformation of the form: 



Pij = -Siirij (25) 



where d is the interplanar distance. The result is: 



Cik = ^Cijkimjmi (26) 

By way of example we may consider the case of crystallographic slip on a {lll}-plane of an fee crystal. 
Then a simple calculation gives d = a/y/3, with a the cubic lattice parameter, and 

Cik = -^~( c n + C44 - ci 2 )Sik (27) 
\Jia 
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(a) Interplanar potential (b) Resolved shear stress 

Figure 1 : Piecewise-quadratic interplanar potential and its derivative. 



In this expression en, C44 and cyi are the three independent cubic elastic moduli. 

If the slip is additionally constrained to take place in the direction of the Burgers vector, as in eq. (|I"5|), 



then the piecewise quadratic interplanar potential fl24[ ) reduces to the form 

<f>(S)=f^\S-m 2 (28) 

where S = S • s and 

C = C ik SiS k (29) 

The function <p(5) and its derivative, which gives the resolved shear stress as a function of slip, are shown in 
Fig-S- 

2.3 Computation of the energy from the phase field 

At zero temperature, the stable configurations of the dislocation ensemble may be identified with the relative 
minima of the energy <|3|> . We begin by attempting a characterization of the absolute mimimizers of (|), and 



defer the more ambitious program of understanding the relative energy minimizers as well. Inserting (17) 
and (28) into (Q) leads to the energy: 



E[(]= inf £[C|£] (30) 
(ex 

where 

C = 8/b (31) 
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is a normalized slip function, and 

fib 2 



mm 



i 



fib 2 



K\(\ 2 d 2 k - I bs(d z x 



In this expression, 



is the resolved shear stress field, and we write 

K 



k>2 



The normalized slip distribution £ of interest now follows from the problem: 

inf MCI 



(32) 



(33) 



(34) 



(35) 



In definition (pOj), the function £ may be regarded as an integer- valued phase field defined on the slip plane 
S, with the property that £(a?) equals the number of dislocations which have passed over the point x G S. 
Thus, the quantized phase field £ describes an ensemble of perfect, or Volterra, dislocations on S. In writing 
([30|) and (35), X x Y denotes the space of slip distributions and integer- valued phase fields of finite energy, 
i. e., the linear space of functions £ : R 2 -> R and £ : R 2 -> Z such that E[£|£] < oo, which is the 
physically relevant space of solutions. 

Next we seek to express the energy of the dislocation ensemble directly in terms of £. To this end, we 
begin by inserting definition (3C) into the variational problem (|35|), with the result: 



inf inf M£|£l 



(36) 



At this point, we may invert the order of minimization in (|36[), which results in the reduced minimum 
problem: 



where 



inf 



inf M£l 



2d K ^ d + (27T) 2 



fib 2 



K\(\ 2 d 2 k 



J 6s£d 2 x| 



(37) 



(38) 



This minimization leads to a linear problem in £ which may be solved analytically by an application of the 
Fourier transform. To this end, we apply Parseval's identity to obtain: 



E[£] = inf 



f 



C&Y { (2vr) 2 
Minimization with respect to £ gives: 



bs*( d 2 k 



(39) 



£ = Co + 



£ 

1 + Kd/2 



(40) 
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where 



Co 



d s 
Jib l + Kd/2 



(41) 



is the slip distribution for £ = 0, i. e., the slip distribution due to the elasticity of the interatomic plane. 
Insertion of ( f40| ) into ([38|) gives the sought energy: 




where 



En 



4 l + Kd/2' 



d 



l + Kd/2 



d 2 k 



(42) 



(2vr) 2 J 2(j, l + Kd/2 



d 2 k 



(43) 



is independent of £ and represents the elastic energy of the interatomic in the absence of dislocations, i. e., 
at £ = 0. Thus, for instance, if s = r = constant, ( |43| ) corresponds to an elastic energy density per unit area 
of (d/2fi)T 2 , which indeed coincides with the elastic energy density of the interatomic plane. Therefore, the 
excess energy 



(44) 



may be identified with the energy proper of the dislocation ensemble. By virtue of (42j), this energy is now 



expressed entirely in terms of the phase field £. If is evident from the form of (42) that, in the absence of an 
applied field, the dislocation energy vanishes for uniform slip, as required by symmetry. Indeed, in this case 
£ = ndoik), where n G Z and #£)(fe) denotes the Dirac delta at k = 0, and, consequently, K£ vanishes 
identically in the sense of distributions. 



As already noted, the phase field £ in (42) takes integer values and, therefore, represents a distribution 



of Volterra or perfect dislocations. It is evident from the form of (42) that the Peierls potential has the 



effect of regularizing linear elasticity so as to render the energy of Volterra dislocations finite. The lattice 



parameter is retained in the energy functional (42) through the interplanar distance d. The effect of the core 
regularization resides in the factor 



<Pd(k) 



1 



l + Kd/2 



(45) 



appearing in both terms of the energy ([42|). This factor tends to suppress wavelengths on the scale of the 
lattice parameter or shorter. The same factor relates the integer-valued phase field £ to the core-regularized 
slip distribution (, eq. (|{]). In this relation, the function tpd{x) acts as a mollifler: the integral of ifdix) over 
the entire plane is one, or, equivalently, ifd{0) = 1; and ifd{x) defines a Dirac-delta sequence, i. e., 



lim ifd(x) = 5 D (x) 

d— >0 



in the sense of distributions. In real space eq. ( |40| ) takes the convolution form: 



C = Co + ¥d * £ 



(46) 



(47) 



which shows that the effect of the piecewise-quadratic Peierls interplanar potential is to smooth out the phase 
field £. A consequence of this smoothing is that the perfect dislocations described by £ acquire a core of 
width commensurate with d. 



9 



M. Koslowski, A. M. Cuitiho and M. Ortiz 



An Exactly Solvable Phase-Field Theory. 



It should be carefully noted that, despite the quadratic appearance of the energy functional (|42|), the 
attendant variational problem is strongly nonlinear owing to the all-important constraint that £ be an integer- 
valued function. Furthermore, the minimization problem is nonconvex owing to the lack of convexity of the 
set X. This nonlinear and nonconvex structure of the variational problem is inherited from the — similarly 
nonlinear and nonconvex — structure of the Peierls potential (^8j). In addition, the variational problem is 
nonlocal, owing to the presence of long-range elastic interactions. These attributes render the problem ( |37| ) 
mathematically non-trivial and confer its solutions a rich structure. Despite these difficulties, the choice of a 
piecewise quadratic Peierls potential lends the problem analytical tractability, as demonstrated subsequently. 

2.4 Illustrative examples 

As an illustration of the type of core structure predicted by the piecewise quadratic model, we may consider 
the case of a general straight dislocation. For simplicity, we restrict our attention to the isotropic case. It is 
convenient to introduce auxiliary orthonormal axes (xi,x' 2 ) with x[ normal to the dislocation line. In this 
coordinate frame one has: 

f (48) 

ik[ 

corresponding to uniform slip over the half -plane x[ < 0. Here, 5d denotes the Dirac-delta distribution, and 
% is the unit imaginary number. A straightforward calculation gives the slip distribution as 

C = 27r ^) (49 ) 

^ ik[(l + ck') y ' 

where 

. 2 _ cos 2 6\ d 

Sm ° + —)-2 (50) 

is a characteristic core width, and 9 is the angle subtended by the normal to the dislocation line and the 
Burgers vector. Thus, the case 9 = corresponds to an edge dislocation, whereas the case 9 = it/ 2 
corresponds to a screw dislocation. A plot of the core profile is shown in Fig. |[ In addition, the energy per 
unit dislocation length follows in the form 

E lib 2 ( . 2n cos 2 9\ , R 

- = ^— sin 2 + - log- (51) 

L 4tt \ 1 — v J c 

where R is an upper cut-off radius. It is clear from this expression that c plays the role of a core cut-off 
radius. The width of the core depends on the orientation of the dislocation line relative to the Burgers vector, 
and it attains its minimum (maximum) value for edge (screw) dislocations. 

This example illustrates how the introduction of a Peierls interplanar potential renders the dislocation 
core structure, including the core width, well-defined. In particular, the high wave number divergence in 
the dislocation energy is eliminated. The specific core structure ( p9| ) is, of course, a result of the assumed 
piecewise quadratic form of the Peierls interplanar potential. It bears emphasis that here the interplanar 
potential is regarded simply as a convenient device for regularizing the equations of elasticity, and no attempt 
is made to model real dislocation core structures. However, it seems reasonable to expect that, for sufficiently 
well-spaced dislocations, the only macroscopically relevant core parameter is the core energy per unit length, 
and that the details of the core structure play a limited role as regards the overall energetics of the dislocation 
ensemble. 
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Figure 2: Core structure of an infinite straight dislocation predicted by the piecewise quadratic model. 



The Peierls potential introduces a lengthscale d into the formulation of the order of the crystal lattice 
parameter. The behavior of dislocation loops may therefore be expected to differ sharply according as to 
whether their size is much larger or, contrariwise, much smaller, than d. These regimes are exhibited in the 
simple example of a circular dislocation loop of radius R. In this case, £ is the characteristic function of the 
circle of radius R, e. g., centered at the origin, and its Fourier transform is 

|= f R e- ikrcos ( e -<t>)rdrde = 2ir-J 1 {kR) (52) 
Jo Jo * 

where J\ is the Bessel function of the first kind, and we have written {x\,x-i) = (r cos#,rsin#) and 



(fci, fo) = (k cos 4>, k sin (j>). Insertion of (52|) into ( p2| ) gives 

E 



^ 2 „ D 2^3, 2 /4fl 2 



d? 



1/2 1/2 
1/2 1 



■1 



(53) 



3 2 

where G 2 ' 4 is the Meijer's G function (25), and we have set v = for simplicity. From the asymptotic 
behavior of the Meijer function it follows that 



E 



R^O 



and 



E 



fib 2 4R 



R 



oo 



(54) 



(55) 



As expected, for fixed d the energy of the loop tends to zero as R — > 0. Conversely, for fixed R, the energy 
diverges as d — > 0, which corresponds to the limiting case of a perfect Volterra loop in a linear elastic 
crystal. Once again, this divergence illustrates the breakdown of linear elasticity in the presence of perfect 
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dislocations and the crucial role played by the core regularization in eliminating that breakdown. The form 



(q5J) for the energy of large loops is consistent with expressions derived from linear elasticity (|24|). However, 



it should be noted that the core width c = d/4 implied by ( |55| ) is predicted by the present theory, which 
stands in contrast to the theory of linear-elastic dislocations wherein the core cut-off radius is an ad-hoc and 
extraneous parameter. 

The issue of metastability and the existence of relative energy minimizers may be illustrated simply by 
the stable configurations of a dislocation dipole under the action of a resolved shear stress r. In order to 
preclude uniform slip distributions, £ may be required to decay to zero at infinity. Suppose that slip occurs 
on the interval —r/2 < x\ < r/2, where r is dipole width. Then, the corresponding phase field is 

£=ltintf (56) 
k\ 2 

For purposes of the present discussion it suffices to consider the asymptotic limit of r S> d. In this regime 
(fl2|) yields the relation 

9E "" 2 r sklkirdki _ bT = n f i_ tT , (57) 



dr 2vr(l - v) J 2vr(l - v) r 

and 



dE /lb 2 1 

-R 7) — j\ n- > r^oo (58) 

oc 2tt(1 — v) c 



From these two expressions we find: 



fib 2 



r 



E ~ E + 7—7- r log - - brr, r oo (59) 

27r(l — v) c 



for some constant Eq independent of r and c. Equilibrium now demands: 



dE nb 2 1 , 

or 2tt(1 — v) r 



The dipole energy at equilibrium is, therefore, 



^b 2 



r 



E~E Q + —p log--l , r^oo (61) 

27r(l — v) \ c J 

where r is the solution of (|60|). The preceding analysis shows that, in the regime under consideration, there 
exist stable dipoles of arbitrary width r, with Eq. (EH) supplying the resolved shear stress which equilibrates 



the dipole. It is evident from ( pl| ) that, for sufficiently large r, the dipole energy is positive and, therefore, 
greater than the dislocation energy corresponding to a phase field £ = 0, which is zero. Hence, the dipole 
configuration is metastable and constitutes a relative minimum of the dislocation energy. 



3 Energy-minimizing phase fields 

The variational problem (p7|) which characterizes the stable equilibrium configurations of the slip plane may 
be regarded as a constrained minimization problem, in which the unknown field £ is constrained to take 
integer values everywhere on the slip plane. This constraint is nonlinear and nonconvex, which renders the 



problem (37) mathematically nontrivial. In addition, the long-range elastic interactions between dislocations 



render the energy nonlocal. These difficulties notwithstanding, in this section we obtain a distinguished set 



of minimizers of the energy (42) analytically in closed form. 
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3.1 Constrained minimizers and energy-norm projection 



A compelling characterization of the solutions of the constrained minmization problem (37) may be derived 
as follows. Let 77 be the unconstrained energy minimizer. Thus, 77 is the solution of the unconstrained 
minimization problem: 



inf Mr/] 



(62) 



In particular, 77 is allowed to take arbitrary real values. Then, the constrained phase field is the solution of 
the problem (e. g., (|26|)): 



inf 

v 2 



where 



V (u, u) 



is an energy norm corresponding to the inner product 



(u,v) 



1 



Lib 2 



K 



2 l + Kd/2 



u* v d k 



The solution of (63) may be expressed as 



(63) 



(64) 



(65) 



(66) 



where Px denotes the closest-point projection of Y onto X in the sense of norm (|64|). 

The solution procedure leading to the core-regularized slip distribution ( may now be broken down into 
three steps. The first step consists of the solution of the unconstrained minimization problem (52), resulting 
in the real-valued unconstrained minimizer 77. The second step consists of the solution of the constrained 
minimization problem (63), resulting in the projection of 7/ onto its closest integer- valued phase-field £. 
The third step consists of the smoothing of £ according to (40), leading to the introduction of smooth cores 
around all dislocation lines. 

The exact evaluation of the projection ( 53 ) entails some difficulty. The corresponding stationarity con- 
dition consists of the requirement that the functional be stable with respect to a unit increment or decrement 
of £ over an arbitrary region of the slip plane. Thus, if R C M 2 is one such region and xr denotes its 
characteristic function, stationarity demands that the inequalities 



Z + XR 
£~XR 



> 
> 



(67) 
(68) 



be simultaneously satisfied. After some straightforward manipulations ( 67 - 53) may be recast in the form 



2 11 XR 



1 

V> Xr) > -5 



XR 



(69) 



These inequalities set bounds on the possible values of the resultant of the residual tractions acting over 
arbitrary regions of the slip plane. 

We proceed to proof the remarkable result that £(x) is the closest integer to rj(x) pointwise on the slip 
plane, i. e., 



£(se) = Pzrj(x), Vcc e 



(70) 
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where Pi denotes the closest-point projection of R onto Z, i. e., Vx G R, P%x is the integer closest to x. In 
order to verify the bounds (|69|), we choose circles of radius e — » as the test regions i?. Let B e denote the 
circle of radius e centered at the origin, and set 



xr{v) = xbXv - x ) 

where x ranges over all of R 2 . Given a function u : R 2 —> R, we have from (65) 

lib 2 K 



(u, Xr) 



1 



u*{k)e ikx X B c {k)d 2 k 



(71) 



(72) 



(2vr) 2 J 2 l + Kd/2 

In order to elucidate the structure of this inner product as e — ► 0, we may introduce the normalized variables: 

x' = —, k' = ek (73) 



whereupon ( 72[ ) becomes 



1 



fib 2 



(2vr) 2 y 2 1 + (K'/e)d/2 



/V "" ^(feOe^-'^Cfc')^ 



In the limit of e — ► we have, asymptotically, 



{u,Xr) 



1 



(27T) 2 

or, in terms of the original set of variables, 



d 



^-u*(k')e ik ' x 'xBAk 



d 2 k' 



In addition, distributionally as e — > we have \B, ~ ne 2 5D, and ([7|) simplifies to 

1 



(u,XR) 



(27t) = 



^u*(k)e ik *ire 2 d 2 k 



This expression finally evaluates to 



(U,XR) ~ ^-TT€ 2 U(X) 



Setting u = £ — r\ we have 



whereas choosing u = with \R as i n (|7lp, yields 



lib 2 



i,2 2 



(74) 



(75) 



(76) 



(77) 



(78) 



(79) 



(80) 



(81) 



in agreement with (|54|). Inserting ( |79[ ) and ( J80D into ( |69[ ) finally gives the stationarity condition 

i > e(as) > ~\, VxGR 2 

whose unique solution is (|70|), g. e. <i. 

This result is remarkable in several notable respects. Firstly, given the nonlocal character of the energy, 
the local character of the projection ( f70| ) is unexpected. Secondly, as we have already noted, problem ([38|) is 
non-convex owing to the non-convexity of the admissible set X. It is therefore remarkable that the problem 
is amenable to an exact, and particularly simple, solution. Finally, it is interesting to note again the crucial 
role played by the core regularization of the energy which, provided that d > 0, results in finite energies 
© and ©. 
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3.2 Unconstrained equilibrium equations 



Next we turn to the unconstrained minimization problem (62). The corresponding Euler-Lagrange equation 
is: 

^Kfj = bs (82) 

It is clear from this equation that all functions of the form fj oc 8d, corresponding to uniform slip distribu- 
tions over the entire slip plane, are solutions of the homogeneous equation. Therefore, it follows from the 
Fredholm alternative theorem that a necessary condition for (|82|) to have solutions is that ,§(0) = 0, i. e., 



s(x)d z x = (83) 
This condition expresses the requirement that the resultant of the resolved shear-stress field be zero. If this 



condition is met, then fj follows from (82) as 



f, = ^-- + 2itC5 D (84) 
Kb n 

where C is an arbitrary constant. An application of the inverse Fourier transform to this expression gives 

v = G*s + C (85) 

where the fundamental solution 



represents an unconstrained slip distribution which decays to zero at infinity and is in equilibrium with a 
resolved shear stress in the form of a Dirac delta applied at the origin. 

3.3 The general solution 

We proceed to collect all the solution steps outlined in the foregoing and to provide an explicit expression 
for the slip distribution £ as a function of the applied field s. Thus, the solution £ follows in three steps, 
namely: s — ► rj — ► £ — > C- The first step yields the unconstrained phase field rj, the second the integer-valued 



phase field £, and the third step returns the core-regularized slip distribution (. Gathering the relations (85), 



(pq) and ( pOp corresponding to each of the three steps just enunciated gives, explicitly, 

t = ip d *Px(G*s) + C (87) 

where the projection Px is given by ( f70| ) and, as noted earlier, we require s to have zero mean, whereupon 
the solution is determined up to an arbitrary integer C*£Z. 

At this point, it is illuminating to revisit the question of uniqueness and metastability, especially in view 
of the existence of relative energy minimizers such as presented in Section [l/i] . In this respect, it should 
be noted that the solution procedure presented in the foregoing returns a unique solution, modulo uniform 
slips. Indeed, the unconstrained equilibrium equations are linear and, within the constraints imposed by 
the Fredholm alternative, determine a unique unconstrained slip distribution rj. The local truncation ( f70| ) 
in turn yields a unique integer-valued phase field £. Finally, smoothing as in ( p7] ) results in a unique core- 
regularized phase field £. It is therefore clear that the solution procedure outlined above selects a 'preferred' 
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slip distribution among a vast array of competitors. Thus, for instance, in regions of the slip plane where the 
dislocation density is small, i. e., where the distance between neighboring dislocations is large compared to 
the core size d, metastable slip distributions may be obtained by inserting small loops in equilibrium with 
the local resolved shear stress within the intervening area between dislocations. Examples of metastable 



solutions of this type have been discussed in Section p.4[ It is intriguing that, in those examples, e. g., 
in the case of a plane slipping under the action of a uniform resolved shear stress, the unique preferred 
solution C = delivers the minimum attainable energy, i. e., it is an absolute energy minimizer. An open 
mathematical question of some interest is whether the preferred solution ( |87| ) does indeed always deliver an 
absolute energy minimizer. 

3.4 Averages and macroscopic variables 

Many of the macroscopic quantities of interest pertaining to the behavior of a slip plane may be recovered 
by taking the appropriate averages. In order to render the operation of taking averages over the slip plane 
well defined, we may simply assume that the phase field is periodic with unit cell 0, i. e., 

({x) = ({x + l 1 a 1 + l 2 a 2 ), (Z\Z 2 )eZ 2 (88) 

for some basis vectors {01,02}. Under these conditions, an application of Orowan's relation gives the 
macroscopic slip strain as 

7 = 7^/ ^ x = Tf<0 (89) 



or, in view of (|40|), alternatively as 



7 = ^Jjd 2 x = 10 (® (90) 

where 

7o = j (91) 

is a reference slip strain. Here I is the slip-plane spacing, which is assumed known, and 70 is the slip strain 
which is attained when one loop sweeps over the entire slip plane. The dislocation line density per unit 
volume is related to the phase field as 

p=](m\) (92) 

Indeed, the integral of | V£| over a region of the slip plane simply measures the total length of all dislocation 
lines contained in that region. Finally, the global equilibrium condition (|3]) may equivalently be expressed 
in the form: 

(s) = (93) 
which requires that the mean resolved shear stress be zero. 
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4 Irreversible processes and kinetics 

The preceding developments have focused on the energetics of a dislocation ensemble in an elastic crystal 
endowed with a piecewise quadratic Peierls interplanar potential. In addition to this energetics, the disloca- 
tion ensemble may undergo inelastic interactions with the lattice, resulting in lattice friction; and short-range 
inelastic interactions with an assortment of obstacles, such as second-phase particles, forest dislocations, 
and others. We assume that these interactions are irreversible and, therefore, kinetic in nature. The essential 
assumption is that the crossing of an obstacle by a dislocation 'costs' a certain energy, regardless of the 
direction of crossing. The energy toll depends on the strength of the interaction and is dissipated, e. g., as 
heat, and irreversibly lost to the system. 



4.1 Variational formulation 

The irreversible dislocation-obstacle interactions may be built into the variational framework developed 
previously by recourse to time discretization (||; ^; ^; f7|; |J). Thus, henceforth we consider a sequence of 
discrete times to, t\, . . . , t n , t n +i, ■ ■ ■ , presume the slip distribution Co at time to to be known, and seek 
to compute the slip distributions C 1 , . . . , Q n , C n+1 , . . . , at all subsequent times. The central problem is to 
determine the slip distribution £ n+1 at time t n+ \ given the solution ( n at time t n and the applied field s n+ \ 
at time t n+ i. To this end, following <Q f|; [?|; [Sj) we introduce the incremental work function: 

W[C n+l\ C] = £ [c n + l] _ E[C] + J f{x) \ C+ l {x) _ C{x)ld 2 x (94) 



where the elastic energy E[Q is given by (p0|), and the second term represents the incremental work of 
dissipation. In this term, the field f(x) > represents the energy cost per unit area associated with the 
passage of one dislocation over the point x, i. e., with a transition of the form ((x) — > ((x) ± 1. Thus, the 
field f(x) represents the distribution of obstacles over the slip plane and is assumed known. By the work and 
energy identity, it follows that W equals the total work supplied to the system during the interval [t n , t n+ i]. 
Eq. (Q) simply states that part of this work is invested in raising the energy of the crystal, whereas the 
remainder of the external work supplied is invested in overcoming the obstacle resistance. The updated slip 
distribution now follows from the minimum principle: 

inf W[C +l \C n ] (95) 



It should be carefully noted that the work function (|94j) depends on the initial conditions ( n for the time 
step, which allows for irreversibility and hysteresis. 

4.2 Min-max formulation 

In order to enable the application of the general solution procedure outlined in the foregoing, we proceed to 
rephrase problem < j95t > as a min-max problem (e. g., (|27|), § 36). We begin by noting the identity: 

Ixl = max Ax, Vi £ R (96) 
A€[-l,l] 

Since f{x) > 0, it follows from this identity that 

f f(x)\C + \x)-C{x)\d 2 x= sup I g n+ \x){C + \x)-e(*))d 2 x (97) 
J \g n+1 \<f J 
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which may be viewed as a statement of the principle of maximum dissipation. The field g n+1 (x), which 
is required to lie within the bounds ±f(x), may be regarded as the reaction-force field exerted by the 



obstacles on the dislocation ensemble at time t n +i. Inserting representation (97) into d95J) and using standard 



properties of saddle points of concave-convex functions (e. g., q27|), Theorem 36.3) leads to the problem 

sup inf { E[C +1 ] - E[C] + [ g n+1 (x)(C +1 (x)-C(x))d 2 x\ (98) 
|o«+i|</ C" +1 6r [ J ' J 



For fixed g n+ (x) the central minimization problem stated in (p8|) is now of the form (35) analyzed earlier, 
and the solution is given by (|87|). The problem which then remains is to optimize the field g n+1 (x) in order 
to maximize the dissipation at the obstacles. 

The optimization of g n+1 (x) may conveniently be reduced to a problem with linear constraints by the 



L[ 5 " +1 ,A ± ] 



(99) 



introduction of Lagrange multipliers (e. g., (|27[), § 28), which leads to the Lagrangian: 

fn mf y {e[C +1 ] - E[C] + J g n+1 (x){C +1 (x) - C{x))d 2 x 
+ J \ + (x){g n+1 (x) - f(x))d 2 x - \-(x){g n+1 (x) + f(x))d 2 x 

The corresponding Kuhn-Tucker optimality conditions are (e. g., (^7|), Theorem 28.3): 

C n+1 (x) - ( n {x) = X+(x) - \~(x) 
g n+1 (x)-f(x)<0, -g n+1 {x)-f{x)<Q 
X+(x) > 0, A-(sc) > 
(g n+1 (x) - f(x))\+(x) = 0, {g n+1 (x) + f(x))\-(x) = 

which must be satisfied simultaneously by the solution. In ( |100[ ), the updated slip distribution £ n+1 follows 



(100) 
(101) 
(102) 
(103) 



from the general solution (87) as 

C n+1 =<p d *P X [G* (s n+1 - g n+1 /b)) + C n+1 
and the global equilibrium condition ( S3 ) now becomes 

J bs n+1 (x)d 2 x = J g n+1 (x)d 2 x 



(104) 



(105) 



which requires that the total applied force be exactly equilibrated by the resultant of the obstacle reactions. 



It is interesting to note how the optimality conditions (|100J-|103|) give mathematical expression to the 
assumed frictional interaction between the phase and obstacle fields, including stick-slip behavior. Thus, the 



equalities (103) prevent A + (a;) and A (a;) from being nonzero simultaneously. If, for instance, Q n+1 (x) — 
( n (x) > 0, then by virtue of ( |l02| ) and ( |100| ) it necessarily follows that A~(a;) = and X + (x) > 0, 
and, by the first of ( [T03b , g n+1 {x) = f(x). Conversely, if ( n+1 (x) - ( n (x) < 0, then ( |i02| ) and (ffOOb 



require that A (a;) > and A + (a;) = 0, and, by the second of (|103|), g n+1 (x) = —f(x). These two cases 



correspond to slip conditions at x, or yielding. If the conditions ( |101| ) are satisfied as strict inequalities, then 
it follows from ( |101[ ) that, necessarily, X + (x) = and A - (a;) = 0, which, in view of ( |100| ), requires that 
Q n+1 (x) — C n ( x ) = 0. This case corresponds to stick conditions, or elastic unloading. The theory is thus 
capable of describing the loading-unloading irreversibility characteristic of plastic materials. 
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4.3 Special case of short-range obstacles 

Next we consider the special case of an obstacle field consisting of a uniform Peierls stress and a distri- 
bution of short-range obstacles. This case merits special attention owing to the fact that the problem can 
conveniently be reduced to the determination of the phase field at the obstacles, which greatly facilitates 
calculations. 

For definiteness, we shall assume the resolved shear-stress field to be of the form 

s(x,t) = T(t) + s (x) (106) 

where r(t) is the applied resolved shear stress at time t and sq(x) is a self-equilibrated field representing the 
long-range elastic stresses induced by the obstacles. For instance, in calculations of forest hardening s (x) 
may describe the action on the slip plane of the long-range stress field of the forest dislocations. In this case, 
in the vicinity of a secondary dislocation sq(x) tends asymptotically to the field of the osculating straight 
dislocation, and, therefore, has the form 

s ~ bfi^Q-, r -> (107) 
r 

for some function A(9) which depends on the geometry of the secondary dislocation. Here (r, 0) are polar 
coordinates centered at the secondary dislocation. The corresponding unconstrained phase field rjQ diverges 
logarithmically, i. e., 

7] ~ B(6)logr, r^O (108) 

for some function B{6). The attendant phase field £o represents an infinite pile up of perfect dislocation 
loops encircling the secondary dislocation. The depth of this pile up is truncated and rendered finite when 
the core of the secondary dislocations is taken into account. Since the exact geometry of the secondary 
dislocations away from the primary slip plane is unknown within the present formulation, one possibility 
that immediately suggests itself is to model so as a random field. The use of random fields to model barriers 
to dislocation slip, or other similar processes such as first-order phase transitions in magnetic systems, has 



been proposed by Hardikar et cd. (28) and Sethna (29) 



For simplicity, we shall additionally assume that the obstacle distribution and, correspondingly, the 
phase field, is periodic with unit cell $7. For the obstacle system under consideration, the obstacle-strength 
field may be taken to be of the form: 



N 



f(x) = br p + fi -Xi), xefl (109) 



1=1 



where r p > is the Peierls stress, X{ and fi > 0, i = 1, ... ,N, are the positions and strengths of 
the obstacles in f2, respectively, and the function tpd( x ) represents the structure of the obstacles. As noted 
earlier, the detailed modeling of the structure of obstacles is beyond the scope of the present work. However, 



the treatment of short-range obstacles as points, corresponding to setting ipd = 5d in (105), inevitably leads 
to logarithmic divergences of the phase field under the obstacles. In order to avoid these divergences the 
obstacles must be endowed with a finite core. Conveniently, however, for distributions of well-separated 
short-range obstacles the precise form of ifi^ plays a role only in the immediate vicinity of the obstacles and 
is otherwise largely irrelevant. We shall simply require that the overlap function ipd * ipj define a Dirac-delta 
sequence as d — > 0, i. e., that 

\im(ip d -k(p d )(x) = 6 D (x) (110) 

d — >0 
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in the sense of distributions. Thus, for points such that | x — Xi \ ^> d the overlap function t/^ * (fd may 
be treated as a Dirac-delta. In addition, in order to avoid divergent integrals we shall require that the prod- 
uct tpd(k)(pcL{k) decay as 1/ A; 2 as A; — > oo. For forest obstacles, a simple choice consistent with these 
requirements is 



(111) 



which simply states that the core structure of the forest obstacles is identical to the core structure of the 
primary dislocations. 

The problem to be solved now follows by inserting (109) into the Lagrangian (99). In order to facilitate 
the solution of this problem we make the ansatz that the obstacle reaction field is of the form: 



N 



g n+1 (x) = < +1 + Y J 9? +l ~ * e « 



(112) 



i=i 



where g™ , i = 0, 1, . . . , N, are constants. This ansatz may be verified a posteriori by checking that all 
optimality conditions ( |100| - 103) are satisfied. Inserting (112) into (p9|) and making use of identities (40), 
(p0|) and ( |110D yields the reduced Lagrangian: 

L[g n+ \ A±] = inf \e[C +1 ] - E[C] + M 9o + \ln + i - In) +E^ +1 (C +1 " O ) 



N 



N 



+ \n\ {A+« +1 - br p ) - A « +1 +br p )} + £ \f(g? +1 -/,)-£ W +1 + /*) 

(113) 



i=l 



i=l 



where the energy functional is given by (|42|). The Kuhn-Tucker optimality conditions (100 103) now 
reduce to: 



and 



7 n+1 - 7 n = A+ - A - 

-5o n+1 - 6r p < 
Aff > 0, Aq > 



g n+i _ bT P < o, 



(9 



9? +1 



ra+1 



er +1 -Q = K-K 

-fi<0, -g? +1 -fi<0 

At > 0, A7 > 
/i)A+ = 0, ( 5 ™+ 1 + /,)A7 = 







with £ = 1, 



iV. 



If 7n+i 7^ 7n> i- e -> m the presence of macroscopic slip, eqs. ( |1 14| - 1117| ) simply require that 

bT P sgn(j n+ i - 7 n ) 



(114) 
(115) 
(116) 
(117) 



(118) 
(119) 
(120) 
(121) 



(122) 



where sgn(x) = x/\x\ is the signum function. In addition, if |<?o +1 | < & rP > then it necessarily follows 
from ( |iT^|lT7l ) that 7n+i = 7n> corresponding to elastic unloading. Thus, the Peierls stress has the effect 
of introducing an initial threshold, or yield point, for plastic activity. 
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In evaluating the optimality conditions (118 - 121), the phase field at the obstacles is computed as 



n+l 



PzV, 



n+l 



(123) 



which is a special case of the pointwise projection (f70|). In addition, the unconstrained phase field follows 
by inserting (112) into (85) and taking periodicity into account, with the result 



7n+1 (*) = -X>r +1 1/ 



G{x - x' + l l ai + l 2 a 2 ) 

(i\i 2 )& 2 



(<Pd * ipd) ix' - Xi)d x > + C n+ 



Specializing this expression at obstacle sites we obtain the linear system of equations 

N 



ji+1 



J2G ijg J +1 + C n+u i = l,...,N 



(124) 



(125) 



which directly relates the unconstrained phase field at the obstacles and the obstacle reactions. In addition, 



the global equilibrium condition ( 105) further reduces to 



(126) 



where we have made use of the identity (so) = 0. Eq. ( |126[ ) simply states that the applied resolved shear 
stress must be exactly equilibrated by lattice friction and the resultant of all obstacle reactions. 



In view of (|110j), the influence coefficients coupling pairs of distinct obstacles in (|125j) may be computed 



as 



^ G(xi - Xj + l l ai + l 2 a 2 ), i^j 

(i 1 ,z 2 )ez 2 



(127) 



In arriving at this expression we have made use of the assumption that the short-range obstacles are well- 
separated, so that they effectively interact as point obstacles. The evaluation of the self-interaction coeffi- 



cients Gu in (127) requires some care, as the obstacles can no longer be treated as points and their assumed 
structure must be taking into consideration. In this case we have 



Gu ~ / G{x)(ipd -k ■tp d )(x)d 2 x + G{l 1 ai + l 2 a2) (no sum in i) 

(zi,p)ez 2 -o 



(128) 



where the introduction of the function ipd describing the obstacle core is essential in order to avoid a diver- 



gent value of Gu. If, by way of example, we choose the obstacle structure ( |lll[ ), then a straightforward 
calculation gives: 



G(x)(tp d -kip d )(x)d 2 x 



' d 



2ir 2 b 



(129) 



In general, this constant is well defined and scales as l/bd provided that, as previously assumed, the product 
ijjd{k)ipd{k) decays as 1/k 2 as k — > oo. It bears emphasis that the sole influence of the choice of obstacle 



structure function if) d on the entire model resides in the constant (|129J). Thus, difference choices of t/j d result 
in different values of the constant, but do not otherwise affect the remainder of the model. 
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4.4 Algorithmic implementation 



The preceding relations provide a complete basis for updating the phase field incrementally. The precise 
algorithm employed in the calculations reported subsequently proceeds as follows: 

1. Assembly and factorization. The calculations are started by computing the matrix Gij using formulae 



( |127p , ( J128D and ( |129[ ). This matrix is symmetric and positive-definite and has dimension N x N. 
The matrix Gij is then factorized. It should be noted that, for a constant distribution of obstacles, the 
matrix Gij remains constant throughout the entire deformation process and its factorization may be 
performed once and for all at the start of the calculations. 



2. Initialization. At time to, the constant Co in (125) is set to zero. In addition, the unconstrained phase 



field r/o induced by the residual tractions sq is set up by superposition of fields of the form ( |108| ), 
suitably truncated at the origin so as to account for the obstacle core. The sign of the superposed 
fields is alternated randomly in such a way that the mean value of sq is zero, as required. 

3. Incremental update. Assume that the state of the slip plane is completely known at time t n , and that 



a new value C n+ \ of the constant in (|125|) is prescribed at time t n+ \. Thus, the constant C is chosen 
as a convenient control or loading parameter. The state of the slip plane is updated by means of the 
following operations: 



(a) Stick predictor. Set the predictor unconstrained phase field r] 
tor reactions: 

N 



rjf, and compute the predic- 



(130) 



i=i 



(b) Reaction projection. Project g™ +1 onto the closest point of the admissible set: \gi\ < fi, i 



1,...,N, according to: 



_n+l 



h, if ~g? +l >h 
if \~g? +1 \<fr 

-f h if g:? +1 <-fi 



(131) 



4. Phase-held evaluation. Once the obstacle reactions g n+1 are known, the unconstrained phase field 



rj n+ i (x) may be computed from (124 ). At the obstacles, the unconstrained phase field follows directly 



as 



N 



V 



n+1 



^ Gij9™ +1 + Cn+1 



(132) 



For x sufficiently distant from all obstacles, we may treat the obstacles as points, whereupon ( |124| ) 
simplifies to: 



N 



Vn+l(x) =J2 



9i + W+l 



(133) 



G(x-x i + l 1 a 1 + l 2 a 2 ) 

i=i |_(/V 2 )ez 2 

Finally, the local projection ([7C|) returns the phase field £ n+ i(x), and < J47| ) gives the core-regularized 
phase field ( n+ i(x). 

5. Macroscopic variables. Finally, the macroscopic resolved shear stress r n+ i and slip strain 7 n +i may 



be computed from ( |126| ) and (|90|), respectively. 
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ii y/y 

(a) Stress vs. strain (b) Dislocation density vs. strain 

Figure 3: Monotonic loading, (a) Applied shear stress vs. macroscopic slip strain, (b) Evolution of disloca- 
tion density with macroscopic slip strain. 

5 Application to the forest hardening mechanism 

In this section, we apply the general framework developed in the foregoing to the forest hardening mech- 
anism. To this end, we consider a slip plane traversed at random locations by fixed secondary or forest 
dislocations and acted upon by an applied resolved shear stress. The objective is to characterize the dislo- 
cation patterns which arise in response to the applied loading, and to determine the effective behavior of 
the system measured, e. g., in terms of the macroscopic slip strain and dislocation line density. In situa- 
tions where the forest hardening mechanism is dominant, the kinetics of the primary dislocation ensemble 
is primarily governed by: the interaction between the dislocations and the applied resolved shear stress; the 
long-range elastic interactions between the primary dislocations; the line-tension effect resulting from the 
core structure of the dislocations; the interaction between the primary dislocations and the long-range elastic 
stress field of the forest dislocations; and the short-range interactions between primary and forest disloca- 
tions such as may result in jogs, the formation of junctions, and other reaction products. Thus, the forest 
hardening mechanism provides a convenient framework for illustrating the range of behaviors predicted by 
the theory. 

5.1 Monotonic Loading 

We begin by considering the base case of a single slip system containing a fixed concentration of point 
obstacles of uniform strength and deforming under the action of a monotonically increasing shear stress. 
We assume a periodicity, and the unit cell £1 is taken to be a square of dimension 1006. We randomly select 
100 points within the periodic cell as the obstacle sites. To these obstacles we assign a uniform strengths 
/ = 10//6 2 . For simplicity, we take the Peierls stress t p = 0, and set Poisson's ratio v = 0.3. 

Fig. H shows the computed stress-strain curve and dislocation density as functions of the macroscopic 
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(a) t/t = 0.00 



(b) r/r = 0.20 



(c) r/r = 0.40 






(d) t/t = 0.60 



(e) t/t =0.80 



(f) t/t = 0.99 



Figure 4: Evolution of the dislocation pattern in response to monotonic loading. Figs, (a)-(f) correspond to 
applied shear stresses r/r = 0.00, 0.20, 0.40, 0.60, 0.80 and 0.99, respectively. 
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(d) r/r = 0.60 (e) r/r = 0.80 (f) t/t = 0.99 

Figure 5: Detailed view of a pair of obstacles as is bypassed by succeeding dislocations. Figs, (a)-(f) 
correspond to applied shear stresses t/tq = 0.00, 0.20, 0.40, 0.60, 0.80 and 0.99, respectively. Individual 
dislocations are color coded for ease of tracking. 
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y/y y/y 

(a) Stress vs. strain (b) Dislocation density vs. strain 

Figure 6: Monotonic loading, effect of variable obstacle strength, (a) Applied shear stress vs. macroscopic 
slip strain, (b) Evolution of dislocation density with macroscopic slip strain. 



slip strain. Fig. ^| shows the corresponding evolution of the phase field £ and the dislocation pattern. In this 
latter plot, the dislocation lines are identified with the lines across which the phase field £ jumps by one, or 
equivalently, with the level contours of 77 + 1/2 at integer heights. The phase field itself counts the number 
of dislocation which have passed over a given point of the slip plane, with proper accounting of the sign of 
traversal. Two well-differentiated regimes, with a transition point at 7/70 ~ 10, are clearly discernible in 
Fig. |3|: a first regime of 'micro-slip' dominated by dislocation bow-out; and a second regime characterized 
by generalized yielding. The saturation value of the resolved shear stress in this latter regime follows from 
global equilibrium as 

1 - 

b ro = bT P + —Y,fi (134) 

' ' i=l 

i. e., the saturation stress is set by the mean obstacle strength. During the micro-slip regime, the value of 
the phase field remains close to zero at the obstacles, which are essentially impenetrable to the dislocations. 
The dislocation loops move reversibly by bowing through the open spaces between obstacles, and by subse- 
quently pinching behind the obstacles, leading to the formation of Orowan loops. As increasing numbers of 
dislocations bypass the obstacles by the bow-out mechanism, the number of Orowan loops surrounding the 
obstacles correspondingly increases, Fig. |j. 

The details of this process are illustrated in Fig. |L which zooms on a particular pair of obstacles and 
shows the evolution of the dislocation/obstacle interactions with increasing stress. In this figure, individual 
dislocations move from left to right and are color coded for ease of tracking. Fig. ||a shows the initial 
configuration. In Fig. |5}), the dislocation has bypassed the obstacles and left Orowan loops in its trail. In 
subsequent frames, new dislocations arrive from the left, bow through the obstacles, and eventually bypass 
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them, leaving behind additional Orowan loops. The ease with which the phase-field representation describes 
these geometrical and topological transitions is quite remarkable. 

At larger applied stresses, the obstacles gradually yield and are overcome by the dislocations. Figs, ^d- 
e correspond to the transitional phase between the micro-slip and the generalized yielding regimes. At 
sufficiently large applied stresses close to the saturation stress to, all obstacles yield and are crossed by the 
dislocations. Under these conditions, the phase field £ increases uniformly without change of shape over 
the entire slip plane. In this regime, the level contours of £ remain ostensibly unchanged, even as their 
height rises steadily, and the dislocation pattern becomes frozen in place, Fig. ||e and |]f. The hardening 
rate correspondingly drops and the stress-strain curve saturates asymptotically, Fig. ||]a. The presence of a 
saturation stage is consistent with observations of fee single crystals undergoing single slip, e. g., during the 
easy glide or stage I or hardening. 

The evolution of the dislocation line density (92) with slip strain is shown in Fig. |3|b. The computed 



dislocation densities are normalized by the reference density 

* = i (135) 

which is the limiting or saturation density corresponding to an arrangement of parallel straight dislocations 
at intervals of b. As may be seen from the figure, the slip plane contains a nonzero dislocation density at 
zero slip strain. This initial density is induced by the long-range elastic-stress field so of the secondary 
dislocations. The process of slip is accompanied by a steady increase in the dislocation line density (p2|). 
Dislocation multiplication and the proliferation of dislocation loops with increasing slip strain are clearly 
evident in Figs. ^|a-f. In the micro-slip regime, the increase in macroscopic slip is accompanied by a steady 
supply of dislocations bowing between the obstacles. This results in an initial parabolic growth rate, p ~ 7 2 , 
Fig. ||b. By contrast, as already noted when r approaches the saturation stress tq the dislocation pattern 
becomes frozen, e. g., Figs. ^ and f, and the dislocation density remains ostensibly constant. This parabolic 
growth regime and the subsequent saturation phase are indeed consistent with observation ([?(]; [3l]). The 
ability of the theory to predict the evolution of the dislocation density, in addition to predicting hardening 
rates, is noteworthy. 

Fig. |6| shows the effect of a variable obstacle strength on the computed stress-strain curve and dislocation 
density evolution during monotonic loading. In order to investigate this effect we consider obstacle strengths 
of the form / = 10/u6 2 (l + uj), with the random variable lo distributed uniformly within the intervals 
[—0.2,0.2] and [—0.5,0.5]. When the obstacle strengths exhibit variability, the weakest obstacles tend to 
yield first and are overcome by dislocations, which are pinned at the stronger obstacles. Consequently, the 
transition between the micro-slip and generalized yielding regimes becomes more gradual. It is interesting 
to note that, since the average obstacle strength is the same in all cases, the saturation stress is not affected 
by the obstacle strength variability. This example illustrates the ability of the theory to account for the 
combined effect of obstacles of different species. 



5.2 Cyclic Loading 

Next we investigate the behavior predicted by the theory under load reversal and cyclic loading. As in the 
preceding example, we assume a periodicity and take the unit cell Q, to be a square of dimension 1006. 
We randomly select 100 points within the periodic cell as the obstacle sites and assign to these obstacles a 
uniform strengths / = Wpb 2 . The Poisson's ratio v is set to 0.3, and, for simplicity, we take the Peierls 
stress t p = 0. We consider the cases of partial unloading, where the system is reloaded after a small 
amount of reverse yielding; and a fully-reversed loading cycle. In both cases, the unloading and reloading 
events bring to the fore the irreversible and dissipative nature of the system, and necessitate full use of the 
incremental theory developed in Section || 
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y/y Q y/y 

(a) Stress (b) Dislocation density 

Figure 7: Cyclic behavior. Labels a-i indicate the loading sequence, (a) Applied resolved shear stress vs. 
average slip (b) Evolution of dislocation density vs. average slip. 

Figs. [7]a and b show the predicted stress-strain curve and the evolution of the dislocation density in the 
cases of fully-reversed loading. The corresponding evolution of the dislocation ensemble is shown in Fig. || 
for the case of fully-reversed loading. The various matching sections of these curves are labeled for ease 
of reference. The interval O-a of the curves in Fig. [7] simply reproduces the monotonic response described 
in the preceding section. Upon load reversal at point a, the system initially unloads elastically, Fig. [7|a, 
segment a-c. In this stage, the applied shear stress continues to push the dislocations against the obstacles in 
the loading direction. However, as the applied load is decreased the force exerted by the dislocations on the 
obstacles correspondingly decreases and, consequently, remains below the obstacle strength. Under these 
conditions, the dislocations remain pinned at the obstacles and simply recoil elastically. This process gives 
the linear unloading branch a-c in the stress-strain curve, Fig. [7|a. 

It should be carefully noted that, owing to the frictional sliding of the dislocations over the obstacles 
that occurs during the loading state, the interaction forces between the dislocations and the obstacles do not 
vanish when the applied shear stress is reduced to zero. Instead, the system is left in a state of permanent or 
residual deformation in which the dislocation/obstacle force system is self-equilibrated. The corresponding 
dislocation pattern for the fully-reversed case is shown in Fig. ^|c. This remanent field is in analogy to the 
self-equilibrated residual stress fields which remain in elastic-plastic solids upon unloading. The combina- 
tion of a reversed applied stress and the residual force field causes the system to yield prematurely in the 
reserve direction, Fig. [7]a, point d, and thus the theory predicts the Bauschinger effect. 

Once reverse yielding commences, the obstacles continue to oppose the motion of the dislocations, 
which now takes place in the reverse loading direction. In particular, the frictional forces at the obstacles 
switch sign relative to the loading phase. This process of switching is clearly apparent in Fig. ^, which 
shows the evolution of the phase field during the loading cycle. Thus, during the loading phase the obstacles 
pull down on the phase field, Fig. ^b-d, causing it to cusp downward at the obstacles, whereas during the 
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(a)r/r = 0.99 



(b) t/t = 0.50 



(c) t/t q = 0.00 







(d) t/t = -0.50 



(e) r/r = -0.99 



(f) r/r = -0.50 






(g) t/t = 0.00 



(h) t/t = -0.50 



(i) t/t = 0.99 



Figure 8: Evolution of the dislocation pattern in response to cyclic loading. Figs, (a)-(i) correspond to t/tq 
= 0.99, 0.50, 0.00, -0.50, -0.99, -0.50, 0.00, 0.50, 0.99, respectively. 
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(g) t/t = -0.30 (h) r/r = -0.60 (i) r/r = -0.99 

Figure 9: Three-dimensional view of the evolution of the phase field during a fully reversed loading cycle, 
showing the switching of the cusps at the obstacles upon unloading and reloading. Figs, (a)-(i) correspond 
to r/r = 0.00 , 0.30, 0.60, 0.99, 0.60, 0.00, -0.30, -0.60, -0.99, respectively. 
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a C 



d f 




Figure 10: Cyclic reloading behavior for various extents of reverse yielding, exhibiting fading memory 
effect. Labels a-i indicate the loading sequence. 



reverse loading phase the obstacles pull up on the phase field, Fig. ^g-i, causing it to cusp upward. 

The stress-strain curve shown in Fig. 10 reveals that the system exhibits 'fading memory'. Thus, when 
the system is reloaded, the stress-strain curve gradually transitions towards the virgin loading curve, and, 
with sufficient reloading, the system eventually 'forgets' the unloading cycle. The extent of this loss of 
memory depends on the extent of reverse yielding. Consider, for instance, the case in which unloading 
is purely elastic and no reverse yielding occurs. This corresponds to reloading, from point b in Fig. [R]. 
Because no reverse sliding over obstacles has occurred at this point, the dislocation/obstacle configuration 
remains undisturbed. Consequently, upon reloading all the obstacles yield simultaneously and the system 
exhibits an abrupt yield point coinciding exactly with the point of unloading, point a in Fig. [R]. Thus, in 
this case no loss of memory occurs and the system exhibits return-point memory in the sense of Sethna et 
al (j29|). This type of behavior should now be compared with that corresponding to a limited amount of 
reverse loading. In this case, the dislocations slip over the obstacles in the reverse loading direction, and 
some of the details of the dislocation/obstacle configuration established at the unloading point d are lost in 
the process. When the system is reloaded, point e in Fig. 10, the obstacles yield anew gradually and not 
all at once, as in the case of a purely elastic unloading, and the system exhibits a gradual transition towards 
the virgin loading curve, Fig. [R], segment e-f. In the case of fully-reversed unloading, Fig. ^, the amount 
of reverse yielding, segment d-e, is large enough that the unloading point is effectively wiped out from the 
memory of the system, with the result that the reloading curve e-i is ostensibly identical to the unloading 
curve a-e. These trends are in good agreement with the experimental cyclic stress-strain data for structural 
steels reported by Ortiz and Popov (32), which was obtained from tests specially designed to exhibit the 
fading memory effect just described. 

The evolution of the dislocation density during a loading cycle is of considerable interest as well, Fig. [7]b . 
Thus, upon unloading the dislocation density decreases as a result of the elastic relaxation of the dislocation 
lines. The dislocation density bottoms out — but does not vanish entirely — upon the removal of the 
applied stress, Fig. ^ b, point c, as some dislocations remain locked in within the system in the residual 
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Figure 1 1 : Evolution of the dislocation density with macroscopic slip strain during cyclic loading for differ- 
ent values of the obstacle density c. 



32 



M. Koslowski, A. M. Cuitiho and M. Ortiz 



An Exactly Solvable Phase-Field Theory. 




(a) Stress (b) Dislocation Density 

Figure 12: Effect of dislocation line-energy anisotropy on the hardening characteristics of the system, (a) 
Applied shear stress vs. macroscopic slip strain, (b) Evolution of dislocation density with macroscopic slip 
strain. 



state, Fig. |§:. The dislocation density increases again during reverse loading, Fig. ^ b, segment c-d, and 
the cycle is repeated during reloading, Fig. [7]b, segment c-d, segment e-i, giving rise to a dislocation density 
vs. slip strain curve in the form of a 'butterfly'. This type of behavior is indeed observed experimentally 
(|33|). It also arises in models of the stored energy of cold work (34), and is in analogy to the hysteretic loops 
exhibited by magnetic systems (35; 29; 36). 

Finally, we exercise the model over the range of periodic cell sizes 10 2 b - 10 5 b, while at the same 
time keeping the number of obstacles constant at 100. The corresponding obstacles density c = iV/|Oj 
consequently varies from 10~ 2 /6 2 to 10 _8 /fe 2 . Of primary interest here is to ascertain how the macroscopic 
behavior of the slip plane depends on the obstacle density. Our calculations show that the stress-strain curve, 
when expressed in terms of the normalized variables t/tq and 7/70 is ostensibly independent of c. It should 
be carefully noted that, for obstacles of uniform strength, to = fc, and thus the saturation strength scales in 
direct proportion to the obstacle density. The variation of the dislocation density evolution with c is shown 



in Fig. H|. It is interesting to note from this figure that the saturation value of the dislocation density scales 
as 1/ y/c, but that otherwise the evolution of the dislocation density with macroscopic slip strain remains 
ostensibly identical in all cases. 

This example also serves to underscore the ability of the theory to effectively deal with large domains 
of analysis. Thus, owing to the absence of a computational grid the complexity of the calculations scales 
with the number of obstacles and is otherwise independent of the cell size, which enables the consideration 
of large cells sizes in the relevant range of observation. 
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(a) 1/ = 0.0 



(b) v = 0.3 



(c) v = 0.5 



Figure 13: Effect of dislocation line-energy anisotropy on the evolution of the dislocation pattern. Figs. 
(a)-(c) correspond to t/tq = 0.99, and 7/70 = 23, 20 and 18, respectively. 



5.3 Dislocation line-energy anisotropy 

We conclude this section with a brief parametric study of the effect of dislocation line-energy anisotropy on 
the hardening characteristics of the system. In particular, we consider systems which are identical in every 



respect except for the value of Poisson's ratio. By virtue of eq. (51), different choices of Poisson's ratio 
result in different ratios between the energy per unit length of edge and screw segments. We recall that, for 
v > 0, the minimum (maximum) energy per unit length is attained for pure screw (edge) segments. All 
remaining parameters are as in the preceding examples. In particular, the unit cell dimension is 1006, the 
number of obstacles is 100, the obstacle strength is / = Wfib 2 , Poisson's ratio v is 0.3, and the Peierls stress 
t p is set to 0. 

Fig. 12 shows the predicted stress-strain curves and dislocation density evolution for v = 0.0, 0.3 and 



0.5. Since the saturation stress to solely depends on obstacle strength and concentration, all three stress- 
strain curves saturate at the same level. By contrast, as v is increased the dislocation line-energy of both 
edge and screw segments increases and the system correspondingly stiffens, with the result that lower values 



of the slip strain are attained, Fig. |12| a. The rate of dislocation multiplication also decreases with increasing 
v, Fig. |l2|b. Perhaps more interestingly, Fig. 13 shows the effect of the anisotropy on the dislocation geom- 
etry. Thus, as Poisson's ratio is increased screw segments become energetically favorable relative to edge 



segments, and the dislocation pattern exhibits a preponderance of elongated screw segments, Fig. 13c. 



6 Summary and Concluding remarks 

We have developed an analytically tractable phase-field theory of dislocation dynamics, strain hardening and 
hysteresis in ductile single crystals at low temperatures. The phase-field representation furnishes a simple 
and effective means of tracking the motion of large numbers of dislocations within discrete slip planes 
through random arrays of point obstacles under the action of an applied shear stress. The theory rests on a 
variational framework for dissipative systems and accounts for energetic and kinetic effects (||; ||; ||; [7]; [8]). 
The energetics accounted for in the theory include: the core energy of the dislocations, represented by a 
piecewise-quadratic Peierls potential (H); the long-range elastic interactions between primary dislocations 
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and between primary and forest dislocations; and the energy of interaction with the applied resolved shear 
stress. The kinetics of the system stem from the assumed irreversible interactions between dislocations 
and obstacles, and from lattice friction, and result in hardening, path dependency, and hysteresis in the 
macroscopic behavior. 

The theory predicts a range of behaviors which are in qualitative agreement with observation. Thus, in 
the base case of single slip under monotonic loading, the theory predicts saturation at a stress which depends 
on obstacle strength and density. It also predicts dislocation multiplication at an initial parabolic rate, during 
a 'micro-slip' regime, followed by saturation at a maximum dislocation density. The ease with which the 
phase-field representation allows for complex geometrical and topological transitions, as the dislocation 
ensemble percolates through the obstacles, is quite remarkable. The cyclic behavior predicted by the theory 
is particularly noteworthy. As required, the theory predicts the essential phenomena of loading/unloading 
irreversibility and hysteresis. But the theory also captures more subtle aspects of the cyclic behavior of 
crystals, such as the Bauschinger effect, consisting of the premature yielding of a slip system under reverse 
loading; the fading memory effect, whereby reverse yielding gradually eliminates the influence of previous 
loading; and the evolution of the dislocation density, leading to characteristic 'butterfly' curves also observed 
in magnetic systems (35; 2S; 36). 

There are other aspects of the theory which help to make contact with actual materials. For instance, 
it is possible within the theory to account for obstacles of different species, and to accord to them varying 
densities and strengths. This furnishes an avenue for building into the theory results from recent work 
concerned with dislocation-dislocation interactions (10; 11; 37; 12). However, in order to build additional 
physical realism into the theory a number of extensions and enhancements immediately suggest themselves, 
to wit: 



1 . It should be possible, within the general framework outlined in Section ||, to formulate a fully three- 
dimensional theory of crystallographic slip accounting for the full complement of slip systems in a 
crystal class and, e. g., an infinite stack of uniformly-spaced slip planes within each slip system. 

2. A topological restriction of the phase-field representation of crystallographic slip, as presented here, is 
that, since the dislocation lines follow the level contours of a scalar phase field, they must necessarily 
trace closed loops or form networks of zero net Burgers vector. This topological constraint precludes 
the effective representation of cross slip, and therefore points to the need for non-trivial mathematical 
extensions of the theory. 

3. In the formulation presented in this paper, the inelastic interaction between dislocations and the crys- 
tal lattice, and, between dislocations and obstacles, has been modeled as a rate-independent frictional 
interaction in the interest of analytical tractability. However, the incremental variational framework 
developed in (^J ^ ^; ^; ||) is general enough to enable consideration of more complex kinetic rela- 
tions, e. g., accounting for the temperature and strain-rate dependence of the Peierls stress. 

4. As presented here, the theory accounts for the edge/screw line-energy anisotropy predicted by linear 
elasticity. However, this type of anisotropy does not suffice to model the vastly different behavior of 
screw and edge segments observed in bcc crystals, which largely owes to differences in the mobility of 
both types of segments. Within the present theory, a proper accounting of this effect can be achieved 
simply by allowing the Peierls stress to depend on the gradient V£ of the phase field, which describes 
the local orientation of the dislocation segments within the slip plane. 



Finally, the fact that the theory permits the reconstruction of the full dislocation ensemble from the 
value of the phase field at the obstacle sites can hardly be overemphasized. Indeed, the geometry of each 
dislocation line is given by the theory analytically and in close form. In addition, the same analytical 
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treatment endows each dislocation with a well-defined core, which effectively regularizes — and eliminates 
the divergences which afflict — linear elastic dislocation theory. The practical effect of this analytical 
tractability is the elimination of any need for the use of numerical grids to discretize the slip plane or, 
worse, the entire crystal. Since the theory targets the equilibrium configurations of the dislocation ensemble 
directly, no costly resolution of transients, often introduced as numerical artifacts in order to advance the 
solution in time and requiring the use of exceedingly small time steps, is required. These attributes render 
the complexity of the calculations required by the theory, in its present state of development, commensurate 
with the number of obstacles. This drastic reduction in complexity may open the way for embedding the 
theory within a finite-deformation formulation of single-crystal elastic -plastic behavior (e. g., (0; [38])), with 
a view to its use in large-scale finite element calculations of macroscopic samples. 
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